Physical signatures of discontinuities of the time-dependent 
exchange-correlation potential 

Daniel Vieira,"'''' K. Capelle"' and C. A. Ullrich^ 

" Departamento de Fisica e Informdtica, 
Instituto de Fisica de Sao Carlos, Universidade de Sao Paulo, 
Caixa Postal 369, Sao Carlos, Sao Paulo 13560-970, Brazil. 
^ Department of Physics and Astronomy, 
University of Missouri, Columbia, Missouri 65211, USA. 
(Dated: February 6, 2009) 

Abstract 

The exact exchange-correlation (XC) potential in time-dependent density-functional theory 
(TDDFT) is known to develop steps and discontinuities upon change of the particle number in 
spatially confined regions or isolated subsystems. We demonstrate that the self-interaction cor- 
rected adiabatic local-density approximation for the XC potential has this property, using the 
example of electron loss of a model quantum well system. We then study the influence of the XC 
potential discontinuity in a real-time simulation of a dissociation process of an asymmetric double 
quantum well system, and show that it dramatically affects the population of the resulting isolated 
single quantum wells. This indicates the importance of a proper account of the discontinuities in 
TDDFT descriptions of ionization, dissociation or charge transfer processes. 
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I. INTRODUCTION 



The most widely used family of exchange-correlation (XC) functionals in density- 
functional theory (DFT) , the local density approximation (LDA) and its gradient- 
corrected semilocal relatives, has been largely responsible for the enormous popularity of 
the theory for electronic-structure calculations. Despite many successes, there are impor- 
tant properties of the exact XC potential that are missed by the LDA and many other 
popular approximations. In this paper, we shall deal with one such property, namely the 
discontinuity of the XC potential upon change of particle number ^. 

It is well known that the exact XC energy functional must vary linearly as a function 
of the total number of particles N, displaying a derivative discontinuity every time the 
system passes through an integral value of p, l6|. The delocalization error 7|], however, 
causes the LDA to predict a nonlinear curvature for the change of XC energy with particle 
number. This has important practical consequences, for instance leading to a description of 
molecular dissociation where the resulting isolated atoms end up with unphysical fractional 
electron numbers. Therefore, attempts to model molecular dissociation processes, or any 
phenomenon that is associated with the transport of charges between well-separated spatial 
regions or subsystems, must take the discontinuity of the XC functional into account 8|]. 

In the framework of time-dependent density-functional theory (TDDFT) [o], |lo| , discon- 
tinuities of the time-dependent XC potential upon change of particle number have recently 
been observed numerically UJ, and shown to be essential for a proper description of 
strong-field ionization processes. In such processes featuring rapid electron escape, the for- 
mation of steps in the XC potential was observed, which eventually develop into fully-fledged 
discontinuities. These previous time-dependent studies dealt with one-dimensional model 
atoms, either putting in the potential discontinuity by hand ll(| and studying its influence on 
single versus double ionization probabilities, or calculating the exact XC potential, including 
steps and discontinuities, by inverting an exact solution of the time-dependent Schrodinger 
equation 

In this paper we shall work along somewhat different lines: our primary interest is in the 
question how these features of the time-dependent XC potential influence the behavior of a 
dissociating system. Rather than carrying out a numerically exact calculation of a simple, 
few-electron model system, we shall work with physically more realistic (yet still simple) 
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systems, namely doped semiconductor quantum wells in effective-mass approximation, and 
study the electron loss of single quantum wells, and the dissociation of an asymmetric double 
quantum well system. In our TDDFT simulations of these processes, we shall compare the 
adiabatic LDA with and without two types of self- interaction correction (SIC). We will first 
show that the LDA plus SIC indeed leads to jumps in the XC potential under the appropriate 
conditions. Then, we will investigate how these jumps affect the dissociation of the double 
well system. 

The paper is organized as follows: Section 2 gives some theoretical background on the 
concept of self-interaction error and its correction in DFT and TDDFT, as well as the basics 
of our quantum well model system. In Sections 3 A and B we present results of time- 
dependent simulations of electron loss and dissociation of single and double quantum wells, 
respectively. Finally, in Section 4 we give our conclusions. 



II. THEORETICAL BACKGROUND 



A. Static and dynamic self-interaction correction 

In static spin-DFT (SDFT), the noninteracting Kohn-Sham (KS) electrons move in a 
spin-dependent effective potential given by 

Va[n]{r) =Ve^t,a{r) + vu[n]{r) + v^cA'^h^lK^)- (2-1) 

The density is obtained via n{r) = X]cr=t j, ^o-l"^) = '^a-=i lYlT"^ ^kai^), where the orbital 
densities are related to the KS orbitals by rikair) = \{pk„{r)\'^. The first term in Eq. fl2.1l) is 
the external potential, followed by the Hartree and the XC potential. 

One of the exact conditions that the XC potential must satisfy is to be free of the 
self-interaction error. This means, in particular, that the Hartree and XC potential must 
cancel each other in a one-electron system. This is automatically true if, for instance, 
we approximate v^c,cr by the exact-exchange (EXX) potential [l^. However, for LDA and 
semilocal functionals this is a hard constraint to be satisfied. Therefore, various procedures 
for removing the self-interaction error of an approximate XC functional have been developed. 
Perdew and Zunger [l^ proposed the most widely known SIC for approximate XC energy 
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functionals: 



occ 



ET''>^. ni] = E:r">^M (Mnu.] + i^xT^^N^., 0]) . (2.2) 

For a one-electron system it is immediately seen that the XC contributions cancel each other, 
and E,^/^^^ becomes equal to the negative of the Hartree energy, thus cancelling i?H- 
The XC potential associated with E^^^^^ , formally defined by 

xrpPZSICl^ „ 1 

<cr"K>^i]w= 17 '\ (2-3) 



is not straightforward to obtain by direct variation with respect to the density, since 



PZSIC 

XC 



is an explicit functional of the orbital densities but only an implicit functional of 13 1. 



Instead, to obtain a local, state-independent PZSIC XC potential one must apply the so- 
called optimized-effective potential (OEP) llSj methodology or one of its simplifications such 
as the Krieger-Li-Iafrate (KLI) approach 16|. The latter is defined as follows: 

^ occ 

^-J'^^^ = ^J^Tm H Iv'fc'.lr)!' {wxc,fc.(r) + [v^il - u^cm] ] + cc, (2.4) 



where 



and the bars denote orbital averages, e.g. "W^^i- = / d^r\ipk(T{j)\^v^J^J [r) . Being self- 
interaction free, the PZSIC formulation also ensures the correct — 1/r asymptotic behavior 
for the XC potential, contrasting with the exponential decay predicted by the LDA. Such a 
feature is particularly important when dealing, e.g., with the delocalization error [tI], since in 
LDA the outermost electrons are not bound strongly enough. Moreover, it has been shown 
that the OEP/KLI formalism is able to predict the expected potential discontinuity when 



the particle number crosses an integer |16 |. 

The PZSIC approach is by no means the only method to get ric 
A more recent version of SIC, proposed by Lundin and Eriksson 



•interaction errors. 



sel 

17|,Q, removes the self- 



interaction directly in the effective potential, which in this approach is given by 

^'xc?f/''K,^d(r) = -VninuMr) + <r°"K - nu„,n,]{v), (2.6) 

where the XC potential acting on the fcth orbital depends on all orbital densities except 
the kth.. Just like PZSIC, the LESIC approach produces a state-dependent XC potential, 



which also requires the implementation of the OEP/KLI procedure in order to obtain a 

local multiplicative KS potential. However, contrary to PZSIC, the LESIC XC approach 

is formulated directly for the potential, instead of for the energy, so that in principle we 

cannot use Eq. (12.51) of the OEP/KLI formalism. Nevertheless, one can assume that there 

exists an energy functional (even though we do not know its structure) which gives rise 

to the LE orbital potential (12. 6p . and adopt 'Uxc,A:cr(r) = v^c,ka'~^[^(^i''^Ai^) ■^'Qs. (12. 4p 

and (12.51) . Other features such as the asymptotic behavior and potential discontinuity are 

also recovered by LESIC, but unlike PZSIC it predicts an erroneous correction even for the 

□ 

hypothetic exact functional [191 ]. 

Let us now consider dynamical situations. In TDDFT, the time-dependent analog of Eq. 
(EH) is 

v„[n]{v,t) = fext,<7(r, t) + t;HM(r,t) +fxc,aK,'^i](r,t), (2.7) 

where n(r, t) = X]o-=t j. XIT'^ lv'fco-(r, ^)P follows from the time- dependent KS orbitals. 
fxc[?^](r, t) is an even more complex object than the static XC potential, since at any given 
time t it must retain memory effects from all previous times t' < t. Modeling this is a 
formidable task, and many commonly used TDDFT approaches therefore make use of the 
adiabatic approximation, which takes an approximate XC potential from static DFT and 
evaluates it at the instantaneous time-dependent density n{r, t) 20|. Accordingly, our previ- 
ous equations (I2.1l) - (l2.6p . dealing with static SIC and KLI calculations, are readily brought 
under the adiabatic time-dependent umbrella by simply replacing (r) — > (r, t). 

Clearly, adiabatic approximations misrepresent the history of the system, because the 
resulting TD functional does not have a memory of the past. One way to build memory in the 
functional has been outlined in reference 2l|. Here, however, we will focus on aspects where 



memory effects are not expected to be crucial. Therefore, all time-dependent calculations 
we present in the following will make use of the adiabatic approach. 



B. Semiconductor quantum wells 

Testing density functionals in a well controlled and modeled environment is of crucial 
importance for the development of (TD)DFT. In this sense, the quasi two-dimensional elec- 
tron gases which occur in doped semiconductor heterostructures are useful since they rep- 
resent an essentially one- dimensional problem: in the effective-mass approximation, only 
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the direction of quantum confinement needs to be analyzed, whereas the lateral or in-plane 
electronic degrees of freedom can be treated analytically. In particular, GaAs/Al^Gai-^iAs 
heterostructures have attracted considerable attention since the effective-mass approxima- 



22 



23 



24|. 



tion for conduction band electrons works very well in these systems 

We assume that such quantum wells have been grown along the z direction and are 
translationally invariant in the x — y plane. Considering only the conduction band, the KS 
eigenfunctions can be written as 

1 



(2.8) 



where A is an area and j is a subband index, and r|| and q|| denote in-plane position and wave 
vectors. The envelope functions fjaiz) satisfy the following one- dimensional KS equation: 

1 



+ Vextiz) + VH[n]{z) + V^cA'^]i^) 



(2.9) 



2m* dz'^ 

with n{z) = X]<T=t i ^'^('^) ~ So-=t J, ~ ^jo-)\fja{z)\'^ /An, and /i being the conduction 

band Fermi level {h = m = 1 from this point on, so that the effective mass m* is a pure 
number; further, the bare electric charge e is replaced by an effective charge e*). The external 
potential Vext{z) is prescribed by the design of the quantum well. To describe the dynamics 
we then propagate the subband envelope functions (Pja{z) using the time-dependent KS 
equation 



1 



+ t^ext {z,t)+ vn [n] {z,t)+ V^c,a N {z, t) 



(2.10) 



2m* dz'^ 

with the initial condition ipj„{z,tn) = ipj„{z). We use LDA+SIC for Vxcaln] {z,t), and the 
corresponding TDKLI potential 



.KLI, 



z,t) 



1 



2n, 



(/i 



25l | can be written as 



C.C. 

(2.11) 

with, e.g., v^cjii^) = / (^^^\Vj^i^^'t)\'^'^xc^J i^^'^)- As we will demonstrate below, the TDKLI 
potential (12.111) retains desirable features of the exact SIC, such as the potential discontinuity 
each time a new subband j is filled and —1/z asymptotic behavior. In the following, the KS 
orbitals are propagated on a real-space grid in real time with a Crank-Nicholson algorithm 
and a time step of 0.05 a.u. 
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FIG. 1: Upper panel: quantum well geometry for studying the effects of barrier suppression. At 
time t = 0, the well depth is abruptly reduced to 1/3 of the initial value. Lower panel: quantum 
well geometry for studying double-well dissociation. The barrier width grows steadily from 0.2 nm 
at t = to 10 nm at t = 250 a.u. 

III. RESULTS AND DISCUSSION 

A. Barrier suppression of a single well 

Mundt and Kiimmel showed for a one-dimensional lithium model atom that the 
time-dependent EXX potential displays humps at its edges as the system loses electrons 
in an ionization process. In the limit of complete ionization, these humps become more 
and more steplike, reaching an asymptotic constant value corresponding to the intrinsic 
derivative discontinuity of the XC potential as the particle number passes through an integer 
value. Here we study an analogous situation, letting electrons escape from a semiconductor 
quantum well. 



7 




FIG. 2: [Color online] Number of electrons in the second subband versus time, illustrating electron 
escape following a suppression of the quantum well depth to one third of its initial value at t = 0. 

We consider a 20nm square GaAs/Alo.sGao.TAs quantum well with conduction band ef- 
fective mass m* = 0.067 and charge e* = e/vTs, and with an electronic density such that 
the two lowest subbands (j = 1,2) are occupied. Initially {t = 0), the difference between 
the conduction band edges of GaAs and Alo.3Gao.7As is taken as 220meV. Once we have 
established the KS ground states ipj^{z,t = 0), for t > we suppress the barrier height 
to one third of the initial value. In other words, we consider a sudden switching process 
which reduces the quantum well depth by 66%, but keeps the square well shape otherwise 
intact (see upper panel of Fig. [1]). Such a process is, of course, highly idealized and would 
be difficult to realize experimentally. In practice, a static electric field could be used to 
suppress the quantum well barrier on one side. 

Applying an absorbing boundary condition, we then let the system evolve such that the 
electrons in the second subband are almost completely ionized, and focus our attention on 
the behavior of the LDA and SIC potentials. We limit ourselves in the following to the 
exchange-only LDA, v^^^^'^[n^,n^]{z,t) = - [G/tt na{z,t)]^^\ applying the TDKLI/PZSIC 
and TDKLI/LESIC approaches. 

Figure [2] shows the number of electrons per unit area in the second subband, N2{t) = 
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(27r)~^(yU — e2a) J \'f2a{z,t)\'^dz, as a function of time. Due to the lowering of the barrier, 
electrons are able to escape, move away from the quantum well and are absorbed at the 
boundary, so that the norm of the envelope functions steadily decreases, / \ip2a{z, t)\'^dz < 1. 
The electrons in the first subband are much more tightly bound, so that the overall ionization 
of the well is almost exclusively a consequence of emptying the second subband. As one 
can observe, for pure LDA the ionization occurs much faster than for LDA+PZSIC and 
LDA+LESIC. This feature is a consequence of the delocalization error [3] and the incorrect 
asymptotic behavior of the LDA: since the electrons are not strongly enough bound to the 
system, they escape more easily. 

Figures [3^, [3)3 and[2t show snapshots of the XC potentials at different times, in the LDA, 
PZSIC and LESIC approaches. Starting with two occupied subbands at t = 0, the second 
subband is almost completely depleted at t = 100 and accordingly the XC SIC potentials 
exhibit a rigid shift between these two situations. As ionization sets in, step structures form 
and quickly migrate away from the quantum well, along with the escaping electrons. Similar 
effects were observed previously in the EXX functional for the one-dimensional lithium atom 



12l | . As expected, the LDA XC potential exhibits neither the formation of the steps nor the 
resulting jump. 

Figure [3] also clearly illustrates the long-range behavior of the LDA and SIC potentials. 
The LDA exhibits the well-known rapid exponential decay; both SIC potentials have a much 
longer spatial range, approaching a — 1/z behavior in the barrier region, sufficiently far away 
from the well. 

Once the appearance of steps and jumps in the time-dependent XC potential has been 
confirmed, the next question is whether and how this affects any physical observables. Pre- 
vious studies [3] were focused mainly on the XC potential itself. In the next section, 
we present an example where the impact of the discontinuities of the time-dependent XC 
potential plays an essential role, namely the dissociation process of a double quantum well. 

B. Dissociation of a double well 

A better understanding and more accurate description of dissociation processes is a central 



issue of much modern research in DFT and TDDFT 261]. Traditionally this problem has been 



approached from a static point of view [5] . Molecular dissociation displays a basic behavior: 
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FIG. 3: [Color online] Snapshots of the time-dependent LDA, PZSIC and LESIC XC potentials for 
the electron escape process shown in Figure [2l The ionization of the second subband produces a 
jump of the XC SIC potentials which builds up via a step structure that moves rapidly away from 
the well. 

the separated atomic systems must have integral charge, and this must be reproduced by any 
accurate XC functional, either in the static or time- dependent case. Unfortunately, most 
local and semilocal functionals fail to satisfy this requirement and lead to final dissociation 
products with fractional charges. 

It has been shown j^, |26| that, as a molecule dissociates, the XC potential develops a 
sharp peak at the bond midpoint followed by the buildup of step structures. However, these 
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insights were obtained in a quasistatic picture, using ground-state calculations in some model 
systems. Here, we use a double quantum well system as a simple time- dependent model of a 
dissociating lieteroatomic molecule. We analyze both the XC potential itself and the effect 
of jumps or steps on the total number of particles placed in each side of a double well system. 

We start at t = with an asymmetric double quantum well divided by an initially very 
thin (0.2 nm) barrier, such that electrons can be viewed as sharing both wells, analogous to 
a molecular orbital. The left well has width 32 nm, and the right one has width 8 nm, and 
both have the same depth of 220 meV (see lower panel of Fig. [1]). The system is populated 
with a given total number of electrons, A^totai = ^left + Alright- 

For t > 0, we let the system dissociate, separating the quantum wells from each other by 
gradually increasing the width of the dividing barrier up to a width of lOnm. The whole 
process is allowed to take a total time of 250 clS indicated in the lower panel of Figure 
[U We chose such a slow dissociation speed in order to avoid the strong density fluctuations 
that would be induced if the wells were torn apart too rapidly; nevertheless, the process is 
far from being adiabatic, and some charge- density oscillations cannot be avoided. At the 
final separation of 10 nm at t = 250 a.u., however, these charge-density oscillations have 
become small ripples, as we will illustrate below. For all practical purposes, the dissociation 
can then be considered complete. 

We monitor the total number of particles in the left and right quantum well, A^ieft and 
Alright, while the system dissociates. In Figure H] we display A^ieft as a function of A^totai when 
the separation of the two wells is complete at t = 250 a.u. The horizontal dotted lines 
indicate the number of electrons A^ieft at which, in a ground-state calculation, the second 
subband in the isolated left quantum well would start to become occupied. This value of 
A^ieft is the same for both SIC calculations, but a smaller value is predicted in LDA. The 
reason for this difference is that SIC (as well as EXX) leads to larger intersubband level 
spacings due to the different asymptotic behavior of the potential |22|] (see also Figure [3]). 

As Figure m shows, the LDA predicts a continuous and rather smooth increase of A^ieft. By 
contrast, the SIC results behave dramatically differently when A^ieft approaches the region 
where the second subband would start to become occupied in the isolated left well. The 
electrons seem to resist filling the left well, and one even sees a marked decrease of A^ieft- If 
A^totai continues to increase, iVieft picks up again, and eventually crosses the threshold of the 
second subband. 
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FIG. 4: [Color online] Number of electrons in the left well, A'loft, as a function of total electron 
number A'totai of a dissociated double quantum well. Both SIC approaches predict an abrupt 
decrease of A'^ieft approaching the region where the second subband would be filled in the isolated 
well (dotted lines). 

In order to understand this behavior, we plot in Figure [5] the time evolution of A^ieft cal- 
culated with PZSIC, up until the dissociated limit of lOnm, with parameters corresponding 
to the points A and B in Figure |H Situation B starts with a slightly larger density in the 
left well, but displays a dramatic drop as soon as the two wells are sufficiently separated. 
The same behavior is found in the LESIC approach, comparing points C and D in Figure HI 
This indicates a repulsive force causing electrons to flow from the left to the right well - a 
clear signature of steps in the potential. 

In Figure [6] we plot snapshots of the PZSIC XC potential and of the density n{z,t) at 
different times, for parameters corresponding to point B in Figures H] and [5l i.e., A'totai = 
3.3 X lO^^cm"^ at the initial time. The snapshots of the XC potential clearly illustrate the 
mechanism preventing the left quantum well from being more filled: Once the double well 
starts to dissociate, the XC potential builds up steps structures, with a very pronounced 
sharp peak in the region between the wells. The system resists putting electrons in the 
second subband: the moment this happens, the potential on the left side shoots up, which 
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FIG. 5: [Color online] Time evolution of Ni^ft for the parameters of point A and B indicated in 
Figure m The drop of particle number in the left well occurs as soon as the system is sufficiently 
separated, so that the subband structure of the isolated left well starts to become relevant. 

means that electrons flow back to the right, as seen from the plots of the density. This effect 
is clearly absent in the LDA. 

It is evident that in a TDDFT computational treatment of molecular dissociation such a 
behavior of XC potential and density must play a crucial role. The key physical requirement 
is the neutrality of the isolated atomic system. By building up steps and jumps in the 
potential, the system avoids incorrect fractional- charge dissociation, forcing electrons to 
flow from an atom to another. This is the main mechanism which restores the nature of 
neutral atoms during the dissociation processes, with the correct integral final charges. 

IV. CONCLUSION 

The results presented here demonstrate that the essential features of the LDA plus SIC 
approach, namely correct asymptotic behavior and discontinuities upon change of particle 
number, appear to carry over from static DFT into the time domain. Comparing two 
different formulations of SIC, PZ and LE, we find only very little difference between them. 
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FIG. 6: [Color online] Snapshots of the PZSIC XC potential corresponding to curve B in Figure [5l 
Due to the rise of the XC potential in the the left well relative to the right, electrons are prevented 
from filling the second subband and are forced to move to the right well instead. A potential 
peak builds up in region between the wells, leading to a repulsive force. The dashed lines indicate 
Vext{z,t), and the dash-dotted lines show the density n(z,t). 



Our TDDFT calculations were performed in the TDKLI approximation, which leads to an 
adiabatic XC potential that is state independent. It is known that the TDKLI approach 
may lead to difficulties with violations of the zero-force theorem 27| , but in our calculations 
this was not a problem (see also Ref. 28l|). 

The time-dependent SIC was used in simulating the electron dynamics in doped semicon- 
ductor quantum well systems. We first studied electron escape of a single square well upon 
sudden reduction of the well depth. The main effect of SIC is here due to its asymptotic 
behavior, which in general makes electron escape proceed at a slower pace compared to 
the LDA. We observed clear indications of a buildup of a step structure in the XC poten- 
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tial, eventually leading to a discontinuity. These observations are in agreement with what 
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Kiimmel and coworkers found earlier for their model atoms 

The most dramatic effect of the XC potential steps and discontinuities was observed 
in the simulated dissociation of an asymmetric double quantum well. This system can be 
compared with a dissociating heteroatomic molecule. Of course, quantum wells are extended 
systems, so that it is not possible to observe the absolute change in particle number. The 
corresponding feature in quantum wells is instead the subband population. We found that 
the dissociation of the double well depends dramatically on the electronic structure of the 
resulting isolated single wells, namely, whether the second subband in the wider well would 
be populated or not. It turns out that in SIC the system tends to resist populating the 
second subband level as long as possible; a pronounced peak structure develops between the 
two wells, and the depth of the XC potential of the wider well jumps with respect to the 
potential in the other well. As a result, the total particle numbers of the left and the right 
well after dissociation are dramatically different from what one would find with the LDA. 

Our calculations were carried out for relatively slowly dissociating quantum wells. Nev- 
ertheless, we were quite far away from the adiabatic limit, as indicated by the presence of 
noticeable fluctuations and charge-density oscillations which, however, subsided once the 
wells were properly separated. The step structures and jumps of the XC potential become 
more and more difficult to discern the more nonadiabatic and abrupt the dynamics becomes, 
and tend to be washed out by strong fluctuations of densities and currents. 

Further systematic studies to identify dynamical regimes and observables for which the 
time-dependent XC potential discontinuities are relevant remain an important task. How- 
ever, there is no doubt that the effect is crucial for generic dissociation or fragmentation 
processes, and can be captured in TDDFT by simple, self-interaction corrected adiabatic 
functionals. 
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